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Sven Nordebo, Senior Member, IEEE, Mats Gustafsson, Member, IEEE, Gerhard Kristensson, Senior 
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Abstract —Natural modes of helical structures are treated 
by using the periodic dyadic Green’s functions in cylindrical 
coordinates. The formulation leads to an infinite system of one¬ 
dimensional integral equations in reciprocal (Fourier) space. Due 
to the twisted structure of the waveguide together with a quasi¬ 
static assumption the set of non-zero coefficients in reciprocal 
space is sparse and the formulation can therefore be used in a 
numerical method based on a truncation of the set of coupled 
integral equations. The periodic dyadic Green’s functions are 
furthermore useful in a simple direct calculation of the quasi¬ 
static fields generated by thin helical wires. 

Index Terms —High-voltage power cables, helical waveguides, 
dispersion relations, open waveguides, volume integral equations. 


1. Introduction 

T he purpose of this paper is to formulate a volume integral 
equation for the determination of the natural modes of 
helical structures. Low-frequency applications are of partic¬ 
ular importance where it can be anticipated the existence of 
twisted modes of a particularly simple structure. The presented 
problem formulation is largely motivated by the need of being 
able to accurately model the field distribution and losses inside 
twisting three-phase high-voltage power cables at 50 Hz, see 
e.g., [8,9]. The approach could also potentially be useful for 
analyzing the wave propagation characteristics of the so called 
litz wires. 

Helical waveguide structures have been treated previously 
such as e.g., with helical sheaths [5,15], and approximations 
for wire helices [16,26]. Presently, there are also very promis¬ 
ing numerical techniques being developed that are based on Fi¬ 
nite Element Modeling [8,9] and the Method of Moments [22]. 
However, to our knowledge there has not been any general pre¬ 
sentation regarding analytical modeling of the natural modes 
of helical structures. It is the aim of this paper to fill in this 
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gap. On the other hand, there is a large body of literature on 
the general dispersion properties of open waveguides, see e.g., 
[12,13,18,19,24], as well as on the general properties of the 
electromagnetic volume integral equations, see e.g., [2^,6, 
25,27]. In particular, it is well known that the volume integral 
operators in electromagnetics are strongly singular and that 
many questions regarding their spectral theory remain largely 
open [6]. Nevertheless, it can be shown in very general settings 
that modes of open waveguides exist and can be interpreted in 
terms of poles of a meromorphic Fourier transform and that 
these poles depend continuously on the model data (except for 
points where poles coalesce or at the boundary of the domain 
of meromorphicity, i.e., at infinity and at the branch-point 
corresponding to the wavenumber of the exterior domain) [19]. 
When there are sources present, the modes (the discrete set of 
eigenfunctions) can be obtained as the residues of the poles 
and the non-discrete set is manifested as an integration along 
the branch-cut [5]. In practical circumstances, the branch-cut 
contribution can often be neglected [17]. 

In this paper, a general helical waveguide structure is treated 
by using classical analytic function theory and Fourier tech¬ 
niques (in particular the convolution theorem) in connection 
with the electromagnetic volume integral formulation and 
a cylindrical vector wave expansion of the related dyadic 
Green’s functions. The Floquet modes are defined as the 
poles of the corresponding integral operators, and the analytic 
periodic Green’s functions are derived by employing the classi¬ 
cal Poisson summation formula. Two independent approaches 
(with and without explicit sources) are used to derive the 
resulting infinite system of one-dimensional integral equations. 
These equations can then be discretized by truncation and by 
using a standard collocation method. 

As a useful byproduct, the periodic dyadic Green’s functions 
provide a simple direct calculation of the quasi-static fields 
generated by thin helical wires. 

H. Natural modes with finite sources 
A. Preliminaries 

Consider a straight helical waveguide structure of radius 
a consisting of a twisting inhomogeneous and anisotropic 
material. The twisting means that the cross-section of the guide 
is rotating along the longitudinal direction of the structure. 
The waveguide may consist of several layers with different 
twist, but it is assumed that the material has a smallest 
common period p in the longitudinal direction. The waveguide 
is assumed to be lossy, and it constitutes an open structure 
placed in a surrounding homogeneous and isotropic free space. 

Let po, eo, po and Co denote the permeability, the permit¬ 
tivity, the wave impedance and the speed of light in vacuum. 
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respectively, and where tjq = \/fJ,ojeo and Cq = l/^Moeo- 
The wavenumber of vacuum is given by kg = ojy/JIoeo where 
w = 277f is the angular frequency and / the frequency. The 
relative permeability and permittivity of the surrounding free 
space are denoted by /i and e, respectively, and the corre¬ 
sponding wavenumber is given by fc = ko^/Jke. The cylindrical 
coordinates are denoted by (p, cj), z), the corresponding unit 
vectors (p, 0, z), the transversal coordinate vector p = pp 
and the radius vector r = p + zz. 

Let E{r) and H{r) denote the electric and magnetic 
fields, respectively, where the time-harmonic factor has 
been suppressed. Further, let Js{r) and Ms{r) denote the 
imposed electric and magnetic sources, respectively, and which 
are assumed to be constrained to a finite region 14 inside 
the waveguide structure. Maxwell’s equations [10] for this 
situation are given by 

V X E{r) = iizpoKr) ■ H{r) - Ms(r), 

V X H{r) = -iu!eoe{r) ■ E{r) -f Js(r), 


finite source region 14, respectively, and where r € 14n, cf., 
[2-6,25]. Here, Ge{r,r',k) and fc) are the electric 

and magnetic dyadic Green’s functions for the surrounding 
free space, respectively, defined by 

Ge(r, r', k) = {I+ ^VV}G(r, r', k), 

G^{r,r',k)=VG{r,r',k)xI, 

ifc I 

and where G{r,r', k) = is the corresponding scalar 

Green’s function, see e.g., [4,5,11]. 

In cylindrical coordinates, the electric dyadic Green’s func¬ 
tion Ge{r,r',k) can be expressed as 

Ge(r, r', k) = G°(r, r', k) - ^ppS{r - r'), (7) 

where G^(r,r',k) is the part of Ge{r,r',k) which can be 
expanded in transverse (solenoidal) cylindrical vector waves 
for r ^ r', d{r — r') the three-dimensional delta distribution 
and ~^pp the corresponding source-point dyadics [4,5]. 
Here, the expansion in cylindrical vector waves is given by 


where /i.(r) and e{r) are the complex valued relative perme¬ 
ability and permittivity dyadics of the material, respectively. 
The Maxwell’s equations (1) are also supplemented with a 
radiation condition providing a unique solution with fields 
vanishing at infinity [4]. The equations (1) can be reformulated 
in terms of the surrounding free space as 

f V X E{r) = iujpok‘H{r) — M{r) — 7Vfs(r), 

{ ( 2 ) 

y V X H{r) = -iu}eoeE{r) + J{r) + Js(r), 

where M{r) and J(r) are the equivalent magnetic and 
electric sources 

I M{r) =-iujpok'Xmir) ■ H{r), 

\ J(r) = -iweoeXe(r) • E{r), 

where Xni(f) and Xeir) are the magnetic and electric suscep¬ 
tibility (or contrast) dyadics, respectively, defined by 

f Xm{r) = lp{r) - I, 

{ (4) 

I Xe(r') = ie(r) - I, 

and where I is the identity dyadic. 

The solution to (1) can now be expressed in terms of the 
following integral equation of the second kind 

Eir) - P Ge(r, r', k) ■ Xeir') ■ E{r') du' 

-iwpoM/y„ Gm{r,r',k) ■ Xm{r') ■ H{r')dv' 

= iujpop fy Ge(r, r', k) ■ Js{r') dv' 

-fy Gm(r,r',k) ■ Ms(r')dv', 

H(r) + iweoe Gm(r, r', k) ■ Xe{r') ■ E{r') du' 

Ge(r,r',k) ■ Xm{r') ■ H{r')dv' 

= fy Gin(r, r', k) ■ Js(r') dv' 

-fiweoe fy Ge(r,r',k) ■ Ms(r')dv', 

(5) 

where the integrals extend over the support of the material 

Vni = {r\{p,(j),z) e [0,a] X [0,27r] x [—oo,-foo]}, and the 


Gl{r,r',k) 

-j pOO ^ 

J — OO ^ — 

( 8 ) 

where the dyadic am{p, p',k,a) is defined in (60) in the 
appendix A. Similarly, the magnetic dyadic Green’s function 
Gni(f,'t'^ fc) is expanded in cylindrical vector waves based 
on the dyadic bm{p,p',k,a) defined in (63). Note that 
Gni(f,'t'^ fc) does not contain any source-point dyadics at 
r = r' [29]. 

To simplify the description below, a non-magnetic material 
is assumed where Xm(f) = 0 and Xe(f) = x(f), and only 
the first equation is needed in (5). A generalization to the full 
system (5) will be straightforward. Extracting the contribution 
from the source-point as defined in (7), the electric field 
integral equation in (5) can now be expressed as 


[1 + PP- X(^)] • E{r) 

-k^ [ GUr,r\k)-x{r')-E{r')dv'= F{r), (9) 

where r G Vm, and the source vector function F{r) is given 
by 


F{r) = iujpoP / G°(r, r', k) ■ Js(r') dv' 

- iujpop^pp ■ Js{r) - J Gmir, r', k) ■ M^ir') du'. 

( 10 ) 

At each radial distance p the material is periodic in the 0- 
coordinate with period 27r and in the z-coordinate with period 
p. The material dyadic x(t’) can hence be represented by the 
following two-dimensional Fourier series 

OO OO 

X{r)= X! X! Xmn(p)e™'^e'”^^ (11) 

m——oo n——oo 

where Xmn{p) are the corresponding dyadic Fourier series 
coefficients. 
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B. Fourier analysis 

Consider the following Fourier representation of the electric 
field 


1 


= ^ / E Ernip,a)e^ 


imd) ^icxz 


da, 


( 12 ) 


' m= —oo 


and similarly for the sources Js{r) and Ms{r) and the source 
vector function F{r). The material dyadic is furthermore 
represented as 


X{r) = 


1 

27r , 


Xm(p,a)e‘ 


da, 


(13) 


' m— — oo 


where 


_^ 27r 

Xm(p, a) =27r V Xnin{p)Sia-n —), (14) 

P 

n— — (X) 

and where the coefficients Xmn{p) have been defined in (11). 

Based on the convolution theorem the Fourier transforma¬ 
tion of (9) yields the following integral equation in Fourier 
space 

^ I'OO 

+ PP'^/ Y Xni-m'{p,a - a) ■ Em'ip,a)da' 
m' — — oo 

Fa ^ Foo ^ 

k 27r / cijji{p-, p , o;)* ^ ^ Xm—m' {p ? ol ) 

® m' — — oo 

■ Em'ip', a') da'p' dp' = Fmip,a), (15) 

where p € [0, a], m G Z and a G M. The Fourier transforma¬ 
tion of the source vector function is given by 

Fmip,a) =\iopqP2ti / amip, p',k,a) ■ J^mip',a)p'dp' 
Jo 

- iw^oPp PP • JsmiP, a) 

-2tt f bmip,p',k,a)-Msmip',a)p'dp'. (16) 
Jo 

By further exploiting the distributional property given in (14), 
the integral equation (15) becomes 

oo oo 

Emip,a)+pp- Y Y Xm-m',liP)-Em'{P,a-l -) 


m ——OO L — — 00 


k 27r / p , k^ ct) • ^ ^ ^ ^ Xm—m',i{p ) 

-'0 / , 

m' — — oo 1——00 

27r 

■ Em'ip',a-l — )p'dp'= F„(p, a). (17) 

P 

The integral operator defined on the left-hand side of (17) 
can be continued to the complex a-plane taking the branch- 
cut (defined by the square root k = s/k'^ — a^) into proper 
account, and hence defines an analytic operator-valued func¬ 
tion Aia) where Aia)E{p,a) = F{p,a). The solution to 
(17) can be expressed as the sum of discrete modes (residues at 
poles), plus an integration along the branch-cut [5]. In practice, 
the branch-cut contribution can often be neglected [17]. 


C. The Floquet theorem 

It is assumed that the sources Js(r) and Ms{r) are sup¬ 
ported only for z < 0, and hence that the source vector 
function Fmip,a) defined in (16) is an analytic function in 
the upper half-plane Ima > 0. According to the Floquet 
theorem [5], there will be Floquet modes in the source-free 
region z > 0, and where the electric field can be expressed as 


OO OO 


Eir)= Y Y Emnip)d''"^e- 


(18) 


m— — oo n— — oo 


where (3 is the complex valued propagation constant of the 
Floquet mode and where Im /? > 0. The Fourier transform 
Emip, a) can hence be written as 


EmiP,^} — ^ ^ Emnip) . 


i(a- (/3-f n^)) 


■^m ip, , 

(19) 

with poles at P+n'^. Here, Emip^ a) is the Fourier transform 
of the left-sided part of Eir) which is supported in 2 ; < 0, and 
hence E^ ip, a) is analytic in the upper half-plane Im a > 0. 
The residue theorem can be stated here as 


1 


/(«)dQ; , Stt 


27r/c„ i(a- ifJ + q^)) 


P 


where C„ is a small closed contour enclosing the point /3+n^, 
/(a) is a function which is analytic inside C„ and Snq denotes 
the Kronecker delta. By applying the residue theorem (20) to 
(17), the following integral equation is obtained for Emnip) 


OO OO 


^mn ip)+pp- Y Y Xm—m',n—n' ip) ■ ^m'n' ip) 

m' — — oo n' — — oo 

27r 

- fc^27r / amip, p',k, (3 + n —) 

Jo P 

OO OO 

^ ^ ^ ^ Xm—m',n—n'iP ) ' Fm'n'iP )P dp = 0, 


m'—— oo n'——oo 


( 21 ) 

where p G [0, a] and (m, n) are integers. 

III. Natural modes with periodic excitation 
A. Preliminaries 

Maxwell’s equations for the free space with a periodically 
modulated plane wave excitation can be written 


V X Eir) = iujpopHir) — Mir), 

V X ff(r) = — iweoe£^(r) 4- J(r), 


( 22 ) 


where the magnetic source is given by Mir) = Af(r)e*^^ 
with a p-periodic part Mir) = Mir -f zp), and similarly 
for the electric source J(r) = J(r)e'^^, and where /3 is 
the Floquet wavenumber with Im/3 > 0. In the transverse 
plane the sources are assumed to be constrained to the cross- 
sectional area S with radius a. The Maxwell’s equations 
(22) are also supplemented with a radiation condition in the 
transverse plane. 
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Following the standard derivation based on vector potentials 
as in e.g., [11], the solution to (22) can be written 

E(r) = iw/xop Jg fg Gep(r, r', k,/3) ■ J (r') dS' dz' 

~ Is lo Gmpir,r',k,^) ■ M{r')dS'dz', 

Hir) = fsfl Gmp(r,r',fc,/3) ■J(r')dS'dz' 

+iweoe fg Gep(r, r', k,/3) ■ M(r')dS' dz', 

(23) 

where the integration is over one unit cell Vc = {r|(p, z) e 
S X [0,p]} and r & Vc, and where Gep{r,r',k, (3) 
and Gynp i'l', r',k, f3) are the periodic electric and magnetic 
Green’s dyadics, respectively, defined by 

r Gep(r,r',fc,/?) = {J+iVV}Gp(r,r',fc,/3), 

\ Gcnpir,r',k,l3)=VGp{r,r',k,f3)xI, 

and where Gp{r,r',k,0) is the scalar periodic Green’s 
function defined in (69) in Appendix B. The peri¬ 
odic electric Green’s dyadic can furthermore be factor¬ 
ized as Gep{r,r',k, P) = Ge{r,r' where 
Ge{r, r', k, /3) is periodic in z—z' with period p, and similarly 
for the periodic magnetic Green’s dyadic Grap{'i',f'',k,0). 


B. Periodic Green’s dyadic 


The spectral representation of the electric Green’s dyadic is 
given by 


Geir,r',k) 


1 . / 

— Ge(p,p',A:,a)e‘“(^"^ Ma, (25) 

J —OO 


where 

Geip,p',k,a) = G°{p,p',k,a) - ^pp6{p- p'), (26) 


and which are based on (57) and (58) in Appendix A. From 
the definition (6) it follows that Ge{r,r',k) can also be 
represented as 


1 roc 1 ^ 

Ge(r, r', k) = ^J {1+^ VV}G(p, p', k, > da, 

(27) 

where G(p, p',k,a) is the scalar two-dimensional Green’s 
function defined in (65) in Appendix B. A comparison of (25) 
and (27) yields immediately that 


Ge(p,p',fc,a)e'“(^-"') 

= {I+ ^VV}G(p, p', k, (28) 

and where both sides can be extended analytically into the 
complex a-plane. Based on the Poisson summation formula 
for the scalar periodic Green’s function given in (74), an 
analytic expression for the periodic electric Green’s dyadic 
can now be derived as follows 

Gep(r, r', k, P) = {I+ ^VV}Gp(r, r', k, /3) 

= ^ fl {I+AvV}G(p,p',fc,/3 + n — 

p p 

n— — oo 

= - £ Gc{p,p',k,P + n—)e'^^+''ik)G-^')_ ( 29 ) 

p ^ p 


Finally, based on (26) the periodic electric Green’s dyadic is 
given by 

Gep(r, r', k, P) = G°p(r, r', k, P) 

- 00 

- ^pp^ip- p') X! Hz - z' (30) 


where 


G%{ry,k,p) 

= - y G«(p,p',fc,/3 + n—(31) 

P P 


Similarly, the following spectral representation of the mag¬ 
netic Green’s dyadic follows from (6) 


1 / 
Gyr,r',k) = — / VG(p,p',fc,a)e‘“(^-^) x/da, 

27r 7_oo 


(32) 


and hence 

G„,(p,p',A:,a)e*“(^-^') = VG(p, p', fc, a)e*“(^-^') x I. 

(33) 

An analytic expression for the periodic magnetic Green’s 
dyadic can then be derived as follows 


Gc^pir, r', k, P) = VGp(r, r', k,p)xl 
= - y VG(p,p',A:,/3 + n —X J 

P P 

n— — oo 

= ^ y Gyp,p',k,p + n—yV+nyv-P)_ ( 34 ) 

p Z—d p 


C. Integral equation for natural modes 

Consider now the source-free Maxwell’s equations as in (1) 
and (2) (with Ms = Jg = 0) and with equivalent sources 
M{r) and J(r) defined as in (3), and where the periodic 
material dyadics are defined as in (11). For simplicity, a non¬ 
magnetic material is assumed here with Xm('f) = 0 and 
Xe(f) = x('t')- The fields are assumed to have a Floquet 
wavenumber P with Im P > 0 and can be expressed as 
E{r) = E{r)e'^^ with a p-periodic part given by 

00 00 

^(^)= E E £^mn(p)e™^e‘”^T (35) 

m——oo n——oo 

The general solution (23) then leads to the following integral 
equation for the electric field 

[T + pp • xir)] ■ E{r) 

-k"^ f f G° {r,r',k,P)-x{r')-E{r')dS'dz' = 0, 

Js Jo 

(36) 

where the integration is over one unit cell 14 and r G Vc, 
and where the contribution from the source-point has been 
extracted by evaluating the periodic distribution in (30). By 



5 


writing G^p{r,r',k,(3) = Gg(r, r', fc ,^ \ the inte¬ 
gral equation (36) can also be written 


[1 + pp- x(t')] • E{r) 

pa p 27 T pp 

-fcM / / G°{r,r',k,/3)-x{r')-E{r')p'dp'd<j)'dz' 

Jo Jo Jo 


= 0, (37) 


The primary interest here is with the twisted modes of 
three-phase power cables at quasi-static (50 Hz) conditions. 
Here, the excitation (16) is governed by a single symmet¬ 
ric component [28] with positive phase progression, i.e., an 
azimuthal Floquet-mode with factor The corresponding 
held twist-mode Emn{p) is dehned by the excitation point 
(to, n) = (1,0) together with the conditions of (21), and which 
yields the following set of feasible Fourier indices 


where r G 14- Here, Gg(r,r',k,i3) is dehned by (31) and 
(59) and can hence be written as 


Gl{r,r',k,P) 


1 °° °° p_ 

= - X! X! a^{p,p',k,l3+n—] 

r> n 






m——oo n——oo 



(38) 


(40) 


where k — 0,±1,±2,... and r = Note that this 

is simply a one-step shift to the right in comparison to the 
illustration in Fig. 1. 


By employing the convolution theorem for two-dimensional 
Fourier series, the integral equation (37) can hnally be written 
in Fourier space exactly as in (21). 


IV. Applications 


A. Twist-modes for multi-conductor power cables 

Assume that the cross-section of the waveguide rotates 
along the longitudinal direction and that there are R radial 
regions TZ^ with distinct rotation (twist) for r = 

In each region, the period in the azimuthal direction is given 
by 271 jmr where to^ G {1, 2,...}. It is furthermore assumed 
that there is a smallest common period p in the longitudinal 
direction such that the longitudinal period pr in each radial 
region is given by pr = pjcir where Ur G {±1,±2,...}, 
and where the positive (negative) sign indicates a left (right) 
handed twist. The twist in each radial region TZ^ is thus 
characterized' by the twist direction {mr^nf), and the non¬ 
zero Fourier coefficients Xmn{p) dehned in (11) can only be 
found at the following points in the reciprocal (Fourier) space 


TO = krUr, 

n = kur, 


(39) 


where fc = 0, ±1, ±2,... and r = 1,..., i?. An example with 
three distinct (one left and two right handed) twist regions is 
illustrated in Fig. 1. 



Fig. 1. Material twist-modes Xmn in Fourier space. Here, T^i, 7^2 and TZs 
denote three twist-regions with (mi,ni) = (3,-1), {m 2 , 712 ) = (7,-3) 
and {m 3 , ns) = (7,3). The filled bullets correspond to a possible truncation. 


'Note that the associated twist angle ‘Pr{p) = arctan(p4^) = 
a.rcta.n{p^ increases with the radius p. 


B. Thin helical wires under quasi-static conditions 

A useful application of the periodic Green’s dyadics ex¬ 
pressed in (23) is with the calculation of the electromagnetic 
helds produced by a helical current distribution. As an exam¬ 
ple, consider the periodic magnetic held produced by a thin 
helical wire carrying the stationary current / under quasi-static 
conditions. The radius of the helix is a and the period is p. 
Here, 

H{r) = I J Gmpir,r',k,0) ■ dr', (41) 

where /3 = 0 and L is a curve constituting one period of the 
helix dehned in cartesian coordinates {x',y',z') as 

( x'= acos{^z'd-ip), 

< (42) 

\ y' = asm{fz' + p), 

where z' £ [0,p] and (p is an off-set parameter of the helix. In 
cylindrical coordinates it is seen that p' = a, f = “^z' -\- p 
and ^ -hi. By using (34), (62) and (63) for p > a, 

the expression (41) yields the result 


Hir) = Ij ^'^fm{p,-rn^)v\^{p',-m^) 

m= —00 r—1 ^ ^ 

■ i^'a— + , (43) 

P 

where the integration over one period has been performed as 


^-i(m+n)^z' 1 / _ 


P ^ dz' = pS-m,n- 


(44) 


V. Numerical method 
A. Discretization by the collocation method 

A discretization based on the collocation method [14] is 
devised as follows. It is assumed that the material region of 
the helical structure is given by the radial domain pi < p < a 
where pi > 0. An A^-point discretization is dehned where 
Pi < P2 <■■■< Pn = o, and Lj (p) denotes the 
corresponding Lagrange basis consisting of linear splines with 
the interpolation property Lj{pi) = 5ij for i,j = 
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see e.g., [14]. The Fourier components of the electric field 
defined in (18) or (35) is now expanded as 

N 

Emnip) = Emno{p)Lj{p), (45) 

J=1 

where Emnjip) is a vector valued coefficient with constant 
cylindrical components. Let pi denote the radial vector p 
evaluated at the point pi fox i = 1,..., N. Due to the 
interpolation property of the Lagrange basis it is seen that 
EmniPi) = Eraniip)- The integral equation (21) evaluated at 
the interpolation points pi now yields the discrete system 


Ejnni ip)+^J2pp'^ m—m' ,n—n' {pi) ' ^m'n'iip) 

m' 7i' 

J a^ipi, p', P + n — ) 

m' n' j — 1 ^ 

\m—m',n—n'{P ')^j{P ) ' E^y^'ji^p )p dp = 0, (46) 

where flj denotes the support region of Lj{p). 

The Fourier series (35) is truncated using M terms based 
on the assumed twist-modes that are at hand, cf., section IV-A. 
A multi-index notation is introduced where k o {m,n), I o 
(to', n') and where k,l = 1,.. ., M. The following definitions 
are made 


Eki Eyyiyiii^p') ^ Elj (47) 

\kil — \m—m'^n—n'^Pi)’, (48) 


and 


f 2 tt 

^kilj (/5) ~ / ^m{,Pi7 P : “f ) 

jQy P 

■ lCm—m',n—n'{_P )Lij{P )P Tp , (49) 


where k,l = 1... ,M and i,j = 1,... ,N. The integration 
in (49) is highly regular with at most some points of dis¬ 
continuity in either of the terms am{pi, p', (3 + n'^) and 
Xm-m',n-n' {p')- It is therefore convenient to place any possi¬ 
ble points of discontinuity of the material function Xmn{p') at 
the grid points pi of the linear interpolation, and to employ an 
efficient quadrature rule based on interior points such as the 
Gauss-Legendre quadrature [14,20] to evaluate the integral 
in (49). The system (46) can now be written in the more 
convenient form 

M M N 

Eki T ^ ^ PP ' Xkil * Eli k ^TT ^ ^ ^ ^kiljiP) ' Elj — 0, 

1=1 1=1 1 = 1 

(50) 

which can be readily interpreted in terms of 3 x 3 block 
matrices and with the row and column multi-indices (fc, i) and 
{l,j), respectively. 

Let A{P) denote the finite matrix corresponding to the 
linear system defined in (50). The propagation constant can 
then be computed from a numerical residue calculus based on 


Pi 


i Ai{A(^)} 


(51) 


where Ai{A(/3)} is the minimum (modulus) eigenvalue of the 
matrix A(/3), and where the closed loop C is circumscribing 
the true zero /3i in such a way that there are no other zeros 
or branch-points of Ai{A(/?)} inside the loop. 

B. Material smoothing 

Appropriate windowing techniques can be applied directly 
in the Fourier domain to reduce the Gibbs phenomena in the 
truncation of the (material) Fourier coefficients Xmn{p)- As an 
example is illustrated in Fig. 2 a comparison of the truncated 
Fourier series expansions of a square pulse with and without 
a one-dimensional coefficient smoothing based on the Kaiser 
window [21]. 


Truncated Fourier series 



Fig. 2. Truncation of the Fourier series of a square pulse, with and without 
spectral smooting and with azimuthal indices m = —20,..., 20. 


Appendix A 

Mode expansions of the dyadic Green’s functions 
Consider a homogeneous and isotropic cylindrical region 
with relative permittivity e, relative permeability p and 
wavenumber k = ki^y/JIe. The solenoidal (source-free) cylin¬ 
drical vector waves are defined here by 

f =-V X , 

< 1 (52) 

[ W2m{r,a) = X wim{r), 

where ipmikip) is a Bessel function or a Hankel function of 
the first kind, each of order to, cf., [1,5]. Here, a is the 
longitudinal wavenumber and k — s/kp — the transver¬ 
sal wavenumber where the square root is chosen such that 
0 < arg K < TT and hence Im k > 0. It can be shown by direct 
calculation that V x a) = kwim{r, a). The following 

curl properties are thus obtained 

W X Wrm{r,a) = kWrm{r,a), (53) 

for r = 1,2, and where f denotes the complement of t (1 = 2 
and 2 = 1). 

The following notation will be used 

Titrm {r,a)= Wrm {p, , (54) 

where r = 1, 2, and the vectors WT-yn{p, a) are given explicitly 
in cylindrical coordinates as 

{ \Tfl 

Wim{p,a) = p -- cjj'lp'yyy^Kp), 

Kp 

io! TTKy. K, 

W 2 m{p, a) = p — 'f'yyy{K,p) - + ZjlljmiKp), 

k ktvp k 

(55) 
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and where the ' denotes differentiation with respect to the 
argument. Let the regular and the outgoing (radiating) cylin¬ 
drical vector waves VTm{r,a) and UTm{r,a) be defined as 
in (52) by using the regular Bessel functions and the Hankel 
functions of the first kind, and {up), respectively. 

The electric dyadic Green’s function defined in (6) can be 
expanded in cylindrical vector waves as 

Ge(r,r',fc) 

/ oo 2 

Y'^rmir<,a)ul^{r>,a)da 

m=-oor=l 

where r< and r> denote the vector in {r,r'} having the 
smallest and largest radial coordinate, respectively, i.e., p< = 
min{/9, p'} and p> = max{p, p'}, cf., [1,4,5]. The dagger 
f refers to a sign-shift in the exponentials in the definition 
(52), and which can be placed on any of the two vector waves 
v^mir, a) and a) in the expression (56), cf., [1]. Note 

that the expression (56) also contains a delta distribution taking 
the source point into account, cf, [5]. 

By employing the notation defined in (54), the electric 
Green’s dyadic can now be expressed as 

Ge(r, r', k) = G°(r, r', k) - ^pp5{r - r') (57) 

where 


and where the dyadic 


bmip,p',k,a) = < 


bm{p,p',k,a) is given by 
J Y '^rm{p, a)vl^{p', a) 

r —1 

J Y '^rm{p, a)u\^{p', a) 

r — 1 


p' < p, 

p' > p- 
(63) 


Appendix B 

Scalar periodic Green’s function 

ifc I 

The scalar Green’s function G{r,r',k) = _^,| satisfies 

the inhomogeneous Helmholtz equation 

{ -I- G(r, r', fc) = —(5(r — r'), (64) 

and can be represented by the Fourier integral 

1 I 

G{r,r',k) = —J G(p, p', fc, Ma, (65) 

where the two-dimensional Green’s function G{p,p',k,a) 
satisfies 


{^1 + K^]G{p,p',k,a) = -5{p- p'), (66) 

and where Vj is the two-dimensional transverse Laplace 
operator, k = \/k'^ — a'^ (ImK > 0) and S{p — p') the 
transverse delta distribution. The function G{p, p', k, a) is 
given by 


G 


1 /-oo 

!(r,r',fc) = — y G°(p,p',fc,a)e'“(^""^ Ma, (58) 


Gl{p,p',k,a) = Y am(p,p',A;,a)e'’”('^ (59) 

m— — oo 

and where the dyadic p' ^k^a) is given by 

. 2 

-Y'^rm{p-,a)vl^{p',a) p' < p, 

T — 1 
. 2 

-Y'"rmip,a)ul^(p',a) p'> p, 


am{p,p', k,a) = { 


T — 1 


(60) 

where the dagger f has been placed on the vec¬ 
tor with primed coordinates, and where = 

The magnetic dyadic Green’s function can be obtained as 
Gni(f, 't', fc) = V X Ge{r, r', k) [11], and it follows from (53) 
and (56) that it can be expanded in cylindrical vector waves 
as 

1 I 

Gm{r,r',k) = — Gm(p, p', fc, Ma, (61) 


where 


Gip,p',k,a) = ^11 ''q\k\P- p'\), (67) 

and which can be expanded as 

. OO 

G{p,p',k,a)='- Y Jr„(NP<)HW(ttp>)e™(^-^'), 

m= —OO 

( 68 ) 

see e.g., [1,11,23] and the Graf’s and Gegenbauer’s addition 
theorem [20]. 

The periodic Green’s function Gp{r,r',k, /3) for plane 
wave excitation satisfies 

{v2 + A:2}Gp(r,r',fc,/3) 

OO 

= —5{p — p') Y^ 5{z — z'+ np)e^^G-^ ) ^ (69) 

n——oo 

together with a radiation condition in the transverse plane, 
and where p is the period and Im /? > 0. The periodic Green’s 
function can furthermore be factorized as 

Gp{r,r',k,P)=G{ry,k,^)e^PG-^'\ (70) 

where G(r, r', k, /3) is periodic in z — z' with period p. 

To derive an analytic expression for the periodic Green’s 
function [23], the following Poisson summation formula [7] 
can be employed 


OO 

Gyp,p',k,a)= Y (62) 

m— — oo 


OO OO 

y] finp)e'^P^ = - Y Hfi + n—), (71) 

^ P ^ P 
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where /3 is real valued and where the Fourier transform is 
defined by 

{ pOO 

F{a)= /(z)e'“^dz, 

(72) 

-| poo ' ' 

= ^ J_^F{a)e-^^^da. 

The periodic samples of G{r,r',k) along the ^'-coordinate 
are given by 

1 ro° 

G{r, r'+znp, k) = — G{p, p\ k, da, 

J —OO 

(73) 

and it follows from the Poisson summation formula (71) that 

OO 

Gp{r,F,k,/3)= ^ G'(r, r' -I-inp, 

n— — oo 

= - y G(p,p',fc,/3 + n— (74) 

p p 

n— — oo 

It is noted that the first sum to the left in (74) converges only 
conditionally (not absolutely) when (3 is real valued and it 
diverges when Im/? 0, whereas the sum on the right-hand 
side of (74) converges exponentially for all (3 with Im /3 > 0. 
By using the property (66) which is valid with complex valued 
wavenumbers k, it is readily seen that the expression in the 
right-hand side of (74) satisfies (69) when [3 is complex valued. 
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